/*!
 ***************************************************************************
 * \file laplace2D.cpp
 * Solving the Laplace equation in 2D.
 ***************************************************************************
 *  \author    Luis M. de la Cruz Salas [ Fri Mar 21 12:26:33 GMT 2010 ]
 ***************************************************************************
 * This code solves the Laplace equation in two dimensions:
 * \f$ \nabla^2 T = 0 \f$.
 * This equation is solved in the unit square \f$ (x,y) \in [0,1]\times [0,4] \f$.
 * The bondary conditions for the temperature are:
 * \f$ T=0 \f$ for \f$ x=0, 4\f$ and \f$ T=10 \sin (\frac{\pi}{2} x) \f$
 * for \f$ y=0, 1\f$
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in
 *  according to your installation and then type the next commands:
 *  \verbatim 
   % make
   % ./laplace2D  \endverbatim
 * 
 * The results are stored in \c Data2D.
 *
 *  The results of this example are shown in the next figures, and are 
 *  plotted using OpenDX with the \c visdat.net script.
 * \image html  laplace2D_01.png "Boundary conditions" width=10cm 
 * \image latex  laplace2D_01.eps "Boundary conditions" width=10cm 
 * \image html  laplace2D_02.png "Final solution" width=10cm 
 * \image latex  laplace2D_02.eps "Final solution" width=10cm 
 ***************************************************************************
 */

#include <iostream>
#include <cmath>
#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Schemes/LaplaceScheme.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/num_utils.hpp"
#include "Utils/inout.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;

// Scalar field in one dimension
typedef TunaArray<float,2>::huge ScalarField2D;

int main()
{
  float length_x = 4.0, length_y = 1.0;
  int num_nodes_x = 40, num_nodes_y = 10; 
  StructuredMesh<Uniform<float, 2> > mesh(length_x, num_nodes_x, length_y, num_nodes_y);
  mesh.print();
  float dx = mesh.getDelta(X);
  float dy = mesh.getDelta(Y);

  ScalarField2D T (mesh.getExtentVolumes());
  ScalarField2D Tn(mesh.getExtentNodes());
  
  firstIndex i;
  Range I(Tn.lbound(firstDim), Tn.ubound(firstDim));
  Tn(I, Tn.lbound(secondDim)) = 10 * sin (0.5 * PI * dx * i);
  Tn(I, Tn.ubound(secondDim)) = 10 * sin (0.5 * PI * dx * i);

  NumUtils::interpolateToVolumes(T, Tn);  

#ifdef WITH_DX
  InOut::writeToFile_DX(Tn, 0, "./Data2D/temp.", dx, dy);
#endif

#ifdef WITH_GNUPLOT
    InOut::writeToFile(Tn, 0, "./Data2D/temp.", dx, dy);
    GNUplot plotter;
    plotter("set grid");
    plotter("splot \"./Data2D/temp.00000\" title \"Initial T \" w l");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    int pause; std::cin >> pause;
#endif
 
    SparseMatrix< Diagonal<float, 2> > A(num_nodes_x, num_nodes_y); 

    ScalarField2D            b(num_nodes_x, num_nodes_y);
    
    ScalarEquation< LaplaceScheme<float, 2> > laplace(T, A, b, 
						      mesh.getDeltas());
    laplace.setDirichlet(TOP_WALL);
    laplace.setDirichlet(BOTTOM_WALL);
    laplace.setDirichlet(LEFT_WALL);
    laplace.setDirichlet(RIGHT_WALL);
    laplace.setGamma(1.0);
    laplace.print();	
    
    laplace.calcCoefficients();
    int iter = Solver::TDMA2DX(laplace, 0.001, 30);
    laplace.update();
    NumUtils::interpolateToNodes(Tn, T);

    std::cout << "\n Iterations = " << iter;
    std::cout << "\t Residual = " << laplace.getResidual();

#ifdef WITH_DX
    InOut::writeToFile_DX(Tn, 1, "./Data2D/temp.", (double)dx, (double)dy);
  std::cout << "\n\n NOTE: The data were written in the [./Data2D] directory in DX format";
#endif

#ifdef WITH_GNUPLOT
  InOut::writeToFile(Tn, 1, "./Data2D/temp.", (double)dx, (double)dy);
    plotter("set grid");
    plotter("splot \"./Data2D/temp.00001\" title \"Final Sol\"w l");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    std::cin >> pause;
#endif

  return 0;
}
